Logarithmic oscillators: ideal Hamiltonian thermostats 



(N 

o 

(N 



O 

B 

I 



I 

O 
o 



> 

00 

\o 

On 

in 

en 

o 

(N 

> 

X 



Michele Campisi, Fei Zhan, Peter Talkner, and Peter Hanggi 
Institute of Physics, University of Augsburg, Universitatsstr. 1, D-86135 Augsburg, Germany 

(Dated: March 1, 2013) 

A logarithmic oscillator (in short, log-oscillator) behaves like an ideal thermostat because ol its 
infinite heat capacity: when it weakly couples to another system, time averages of the system 
observables agree with ensemble averages from a Gibbs distribution with a temperature T that is 
given by the strength of the logarithmic potential. The resulting equations of motion are Hamiltonian 
and may be implemented not only in a computer but also with real- world experiments, e.g., with 
cold atoms. 

PACS numbers: 02.70.Ns, 05.40.-a 67.85.-d 



Thermostats play an important role in computational 
physics [I] . They provide effective and useful methods to 
simulate the action of a thermal environment on systems 
of physical and chemical interest. Mathematically speak- 
ing, their salient feature is to produce "thermostated dy- 
namics" of the system of interest: that is, they are meant 
to impose long-time averages of system observables that 
coincide with Gibbs-ensemble averages at a given tem- 
perature T. Widely used thermostats are: the Langevin 
thermostat [2], Andersen's stochastic collision thermostat 
[3] and the Nose-Hoover deterministic thermostat |3HS]- 

Here we present a thermostat differing in various re- 
spects from the previously reported ones. Our main re- 
sult is that a logarithmic oscillator (or a "log-oscillator" 
as we shall call it below), weakly coupled to the system 
of interest (in short "the system" in what follows) leads 
to thermostated system dynamics. In its simplest ID 
version the system-|- log-oscillator Hamiltonian reads: 
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where p — (pi, . . . ,pn), Q — (<Zi, ■ • ■ ,Pn), fni are the mo- 
menta, positions, and masses of the particles composing 
the system; X, P, M are the log-oscillator position, mo- 
mentum, and mass, respectively; 5 > sets the length 
scale of the log-oscillator and T is the thermostat tem- 
perature; V{q) is the system inter-particle potential and 
h{q, X) denotes a weak interaction energy that couples 
the log-oscillator to the system. When the total Hamil- 
tonian H is ergodic, the system-|-log-oscillator trajectory 
samples the microcanonical ensemble, and the system 
trajectory samples the canonical ensemble at tempera- 
ture T. This continues to hold if the ID log-oscillator 
is replaced by higher dimensional log-oscillators, for ex- 
ample for a charged particle in the attractive logarithmic 
2D Coulomb field generated by a long charged wire. 

Compared to the previously reported thermostats the 
present thermostat exhibits an evident advantage. The 
Hamiltonian (n]) or its higher dimensional versions can be 
readily implemented in a physical experiment. In Fig. [l] 
we show a possible implementation. The system is com- 
posed of a gas of neutral atoms confined into a box. The 



thermostat is an ion subject to the attractive 2D coulomb 
potential generated by a thin oppositely charged wire, 
\QX\/2TTeo In p. Here Q > is the charge of the ion, A < 
the linear charge density of the wire, p the distance be- 
tween wire and particle, and Eq the electric permittivity 
of vacuum. Through short-range repulsive interactions 
the ion thermalizes the neutral gas to the temperature 
T — \QX\/7TeQ. Another possibility for the realization of 
a log-oscillator is by means of a laser beam with an inten- 
sity profile of logarithmic form coupled non-resonantly to 
an atom [7]. This could be realized to thermostat cold 
atomic gases [S]. 

Atomic systems in isolation from the environment nat- 
urally sample the microcanonical ensemble. For small 
systems this sampling may considerably differ from the 
canonical one and can result in distinctive thermody- 
namic features such as negative specific heats. These 
were experimentally investigated with small atomic clus- 
ters [HI [13 • Typically it is difficult to have a small iso- 
lated system sample the canonical Gibbs distribution. 
Our method opens this possibility. More generally, by 
using a single log-oscillator as an environment simulator, 
our method allows to experimentally study thermostated 




FIG. 1: (Color online) A single ion (black sphere) carrying the 
charge Q > 0, subject to the attractive 2D Coulomb potential 
generated by a wire (green cylinder) carrying the linear charge 
density A < 0, thermalizes, by means of short-range repulsive 
collisions, a gas of neutral atoms (orange spheres) to the Gibbs 
distribution of temperature T = {QXI/ttsq. 



small systems in isolation from the real environment. 
One advantage of this paradigm is that our method would 
allow to control a thermal parameter, the temperature T, 
by means of mechanical parameters, e.g., with reference 
to Fig. [T] the charge density A on the wire. 

Just like the Nose-Hoover thermostat, our thermostat 
is deterministic and time-reversible, but at variance with 
Nose-Hoover dynamics which are not Hamiltonian [Tllllj. 
our thermostated dynamics are manifestly Hamiltonian. 
There exist "generalized Hamiltonian formalisms" [1] 
for the Nose-Hoover dynamics in the literature. The 
most prominent examples use Nose's Hamiltonian [1]: 
H^ose - EP?/2m.X2 + V{q) + P^ /2M + fTlnX or 
Dettmann's Hamiltonian [TSJ [T^: Hjj ~ Xi^Nosc- At 
variance with our Hamiltonian in Eq. (fl]), these involve 
the non standard kinetic terms, pf/2miX^ andpf/2miX, 
respectively, which, due to the dependence on the log- 
oscillator position, cannot readily be realized in an exper- 
iment. Further, while Dettmann's Hamiltonian produces 
thermostated trajectories only for a specific value of the 
energy (i.e., Hu = 0) our method thermalizcs the system 
irrespective of the energy value. We elucidated these is- 
sues further in Ref. [13]. The usefulness and importance 
of the Nose-Hoover equations as a computational ther- 
mostat are beyond question ITS]. 

Theory.- Before we shall provide the formal argument 
we present a physical explanation indicating why it is 
plausible that the Hamiltonian in Eq. (fTl) leads to ther- 
mostated system dynamics. Consider the isolated ID 
log-oscillator: 
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Applying the virial theorem, (PdHiog/dP) = 
{XdH\og/dX), to the ID log-oscillator, wc obtain 
(^P^ /M') = T where (•) denotes the time average. That 
means that all trajectories of a log-oscillator have the 
same average kinetic energy [^, i.e., the same kinetic 
temperature (P^/Af) = T, regardless of their energy 
E. This implies dT/dE = 0. Recalling the definition 
of heat capacity, C = dE/dT, one finds that the 
log-oscillator exhibits a spectacular property: its heat 
capacity is infinite, which is the defining feature of an 
ideal thermostat. Since the log-oscillator may only exist 
in the state of temperature T, we expect that a system 
will reach this same temperature T when it is weakly 
coupled to the log-oscillator. 

To formally prove that the log-oscillator induces ther- 
mostated dynamics of the system at the temperature T, 
we recall the general expression for the probability den- 
sity function p(q, p) to find a system at the point (q, p) 
of its phase space when it is weakly coupled to a second 
system [the log-oscillator in the present case], provided 
that the compound system probability distribution is mi- 



crocanonical. It reads [TB] 

niog[Etot- Hs{q,p)] 
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where E'tot is the total (conserved) energy of the com- 
pound system. With E denoting the log-oscillator en- 
ergy, 

^iog(S) = / dXdP6[E - iJiog(X, P)] (4) 

is the density of states of the log-oscillator, and 

^{Etot) = / dXdPdqdp 6[Etot - Hiq, p, X, P)] (5) 

is the density of states of the compound system. Here 
6{...) denotes Dirac's delta function and Hs is the system 
Hamiltonian. 

According to Eq. (|3| the density of states of the log- 
oscillator defines the shape of the distribution of the sys- 
tem. Performing the integration in Eq. Q with the log- 
oscillator Hamiltonian, Eq. ^, one obtains for the den- 
sity of states of the log-oscillator the expression 



ftiogiE) ^ 26V27rA//re^/^ . 
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Inserting Eq. ^ into Eq. ((sl) yields the Gibbs distribu- 
tion for the system, 



p{ci,p) = e-"-^^'P^/^/Z{T). 



(7) 



regardless of the energy Etot assigned to the compound 
system. Here Z(T) = / dqdpe-^-s(q,p)/T jg ^^le system 
canonical partition function. 

Also a /-dimensional log-oscillator _ff(X,P) — 
PV(2A/) -I- /r/21n(XV62) [.^^^ere X and P are vec- 
tors of size /] results in the exponential density of states 
fl\og oc e^'^. Therefore, /-dimensional log-oscillators in- 
duce thermostated dynamics as well. 

So far we have left the system-thermostat interaction 
/i(q, X) unspecified. As in standard statistical mechanics 
where a heat bath with many degrees of freedom replaces 
the single log-oscillator [16], h{c[,X) must comply with 
two requirements: (i) it must be sufficiently weak that 
it can completely be neglected in the calculation of the 
probability density p(q, p). This assumption guarantees 
the applicability of Eq. (l3| provided that the total system 
stays in microcanonical equilibrium. In order that this 
equilibrium state actually is reached from arbitrary ini- 
tial conditions it is necessary (ii) that the total dynamics 
is ergodic. To meet these two requirements, short-range 
repulsive interactions typically suffice, see the numeri- 
cal examples below. Note that with a short-range repul- 
sive interaction, the fraction of time during which the 
log-oscillator interacts with any other particle is much 
smaller than one. This assures that the average inter- 
action energy represents only a small part of the total 



energy, and hence the weak couphng assumption imphed 
by Eq. ^ is met. 

Numerics.- In order to corroborate our statement we 
performed ID and 3D molecular dynamics simulations 
using symplectic integrators [Ifj . 

In our first numerical experiment we used two point 
particles of mass to in a ID box of length L and placed a 
log-oscillator of mass M and strength T between them, 
see the inset in Fig. [2j The three particles interact with 
each other and with the fixed walls via the truncated 
Lennard- Jones potential, reading 



VLj{q) = 
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(8) 

that is h{qi,q2,X) = J^t'^LjilQi - X\), and V{qi,q2) = 
VLj{\qi-q2\)+E,{VLj{h+L/2\)+VLj{h~L/2\)]where 
L is the box length. In the simulations we adopted m, a 
and £, as the units of mass, length, and energy, respec- 
tively. In order to avoid the singularity of the logarithmic 
potential at the origin we replaced it with the following 
potential: 



MX) 






(9) 



For all simulations we used the value h = a. This trun- 
cation results in a correction of the density of states (|6]), 
which vanishes as the energy ii^tot increases. Fig. [2] de- 
picts the probability density function, q{Es) of finding 
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FIG. 2: (Color online) Probability density function of energy 
for a system of two particles in a ID box performing short 
ranged collisions with a log-oscillator. The system energy Es 
is rescaled by the total simulation energy which is Etot ~ 75e. 
The log-oscillator strength is T = 15e and the box length is 



10 e 



Etot/T^ 



1484(7. Black dots: numerical simulation. 



Red line: Gibbs distribution at temperature T = 15e. The 
total system is schematically represented in the inset with the 
system of interest (two orange particles) confined to the box 
potential (orange curve), and the log-oscillator (black parti- 
cle) confined to the logarithmic potential (black curve). 



the system consisting of the two orange particles depicted 
in the inset at the kinetic energy Es in a molecular dy- 
namics simulation at total energy Etot ■ According to Eq. 
^ this should be of the form g{Es) oc e'^^^^nsiEs) oc 
g-Es/T ^ where ^s{Es) is the system density of states. 
Note that ils{Es) is constant in the case of a system 
Hamiltonian composed of two quadratic degrees of free- 
dom. The numerically computed curve excellently fits 
the desired canonical distribution with the expected tem- 
perature T. The simulation energy .Etot was chosen large 
enough, so that the error introduced by the replacement 
of the purely logarithmic potential with the truncated 
one, was negligible. The box length was taken such that 
it exceeded the maximal excursion of the log-oscillator 
a^max — 2o-\/e^^t°t/'^ — 1. Otherwise the log-potential 
would be effectively cut-off by the box-potential and con- 
sequently the exponential shape of the density of states 
would be destroyed. 

Our second numerical experiment considers as ther- 
mal bath a charged particle in the electric field generated 
by a long and oppositely charged wire: the so-called 2D 
Coulomb potential which is of logarithmic form. Fig. IT] 
The charged particle Hamiltonian reads: 
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(10) 
where T — \Q\\/'Ket). Assuming that the motion is con- 
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FIG. 3: (Color online) Probability density function of the 
absolute value of the velocity components of a gas of 3 neu- 
tral particles in a confining box, weakly interacting wi th a 
charged particle in a truncated 2D Coulomb field, Eq. (111. 
The velocity is rescaled by the maximum possible velocity 
Umax = {2Etot/rnY'^ . The simulation energy is Etot = 120e, 
T = 15e, b = a, and the box has dimensions L^ — 20a, L^ = 
Ly = 2rmax (r-max = o-(e^'°''"^ - 1)^/^ ~ 109(T is the maximal 
possible excursion of the log-oscillator). Black dots: numer- 
ical simulation. Black dashed line: Maxwell distribution at 
temperature T = 15e. Red solid line: corrected distribution 
accounting for the truncation of the log-potential Eq. (111. 
Inset: same, but for a gas of 8 particles. 



fined in the Z direction by two rigid walls parallel to 
the XY plane and separated by a distance Lx, one ob- 
tains for the density of states the expression fliog(-E) = 
^5/2(2iVf)3/2L^62rV2e'E/T_ Thus we expect the system 

to behave as a thermostat. In our simulation we let this 
thermostat weakly interact with a neutral gas of 3 parti- 
cles confined in a box, and recorded the probability v(v) 
to find the absolute value of any of the 3x3 velocity 
components of the neutral gas at value v during the sim- 
ulation. As with the ID simulation, the 3-1-1 particles 
were interacting with each other and with the fixed box 
walls via the truncated Lennard- Jones potential, Eq. (|8]) . 
The logarithmic potential is truncated in the same way 
as in the ID case, Eq. (|9]), that is we used the potential 



(^(X, Y) = Tln[(X2 + r^ + \?)l\? 



(11) 



The results are displayed in Fig. [3j The truncation of the 
logarithmic potential entails a deviation of the density of 
states from the exponential form: r2iog(£') oc e^l^\-^ — 
2r(3/2, i?/T)] where r(a,x) is the upper incomplete 
gamma function. Note that with EjT ^ 1 this devi- 
ation vanishes exponentially as Q.\o^[E) oc e^'^ {\/t: — 
2^ E jTe^^l^)^ where we have used the asymptotic ex- 
pansion of the upper incomplete Gamma function J18j . 
This leads to a deviation of the distribution p(w) from 
the Maxwellian form. For a fixed simulation energy E^tot , 
this deviation in p(w) becomes more pronounced as the 
number of degrees of freedom composing the system in- 
creases, cf. the inset in Fig. |3] This can be compensated 
by increasing the simulation energy -Etot- We estimate 
that this scales as Etot ^ c'iNT/2 ~ c{Hs), with some 
constant c depending on the required degree of approxi- 
mation. 

Remarks. - Not only can logarithmic potentials be gen- 
erated artificially, e.g., with properly engineered laser 
fields [7], electrophoretic traps [TH] or with charged wires, 
but they also occur naturally in various situations: For 
example logarithmic potentials govern the motion of stars 
in elliptic galaxies |20j , determine the interaction of vor- 
tices in flow fields [21) . and of probe particles in driven 
fluids [22] • Log-oscillators recently received much atten- 
tion in regard to their anomalous diffusion properties |23l - 
m] . The present work is complementary to these studies 
[23l - f26| in the sense that our focus is on the dynamics of 
the particles surrounding the log-oscillator, whereas their 
focus is on the dynamics of the log-oscillator itself. 

One of the earliest thermostats was proposed by An- 
dersen [3| . In the method of Andersen the system evolves 
according to Hamiltonian equations of motion until, at 
some random time r, the velocity of a randomly chosen 
particle in the system is instantaneously assigned a new 
value drawn from a Maxwell distribution with the desired 
temperature. The system then continues it Hamiltonian 
motion until the next random event occurs, and so on. 
Our method can be seen as a fully deterministic version of 



Andersen thermostat, where the times at which the col- 
lisions occur and the newly imparted velocities are not 
drawn randomly, but follow deterministically from the 
total system dynamics. 

In many studies thermal baths are modeled as infinite 
collections of harmonic oscillators or free particles. In 
the present method this infinite collection is replaced by 
a single log-oscillator. It has therefore the evident ad- 
vantage of not involving any thermodynamic limit while 
retaining the Hamiltonian structure. Roughly speaking, 
the thermodynamic limit is lumped in the singularity of 
the log-potential. At variance with infinite thermal baths 
whose temperature is given by the bath's energy per de- 
gree of freedom, log-oscillator thermostats contain the 
temperature as a parameter in the total Hamiltonian. 
This opens the possibility, for example, to study the re- 
sponse of a system to a varying temperature, and take 
advantage from the non-equilibrium statistical mechani- 
cal machinery dealing with time dependent Hamiltonians 

Another advantage of our method is that, because 
the Hamiltonian is written in the standard physical sys- 
tem+bath+intcraction form: H = Hs + Hb + h, it pro- 
vides a direct way to control the strength of the interac- 
tion h, allowing also to simulate thermalization to gener- 
alized Gibbs states occurring when the system-bath cou- 
pling is not weak [5H], which can be a relevant case for 
small systems. 

Conclusions.- We demonstrated that log-oscillators 
possess infinite heat capacity, i.e., they are ideal ther- 
mal baths. As such they have a thermostating influence 
on the dynamics of many-particle systems. The result- 
ing deterministic Hamiltonian dynamics are distinct from 
the Nose-Hoover dynamics. Unlike previously reported 
generalized Hamiltonian formulations of Nose-Hoover dy- 
namics, our Hamiltonian (i) produces thcrmostated dy- 
namics irrespective of the energy value and (ii) presents 
the kinetic terms in standard form. Consequently it is 
amenable to experimental realization. Its most promis- 
ing practical use is as an analog thermostat simulator for 
the experimental investigation of the thermodynamics of 
small systems, e.g., atomic clusters. 
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